#include <iomanip>
#include <boost/math/distributions/chi_squared.hpp>

#include "cds_vom.h"

#define VOM_BEGIN namespace vom {
#define VOM_END }

VOM_BEGIN

node::node(const str_t* alphabet, int periodic, char ch, node* parent)
  : alphabet_(alphabet), periodic_(periodic), ch_(ch), 
    parent_(parent), children_(alphabet_size(), NULL), 
    count_(periodic, scalar_vec_t(alphabet_size(), 0)), 
    count1_(periodic, scalar_vec_t(alphabet_size(), 0)), 
    f_(periodic, scalar_vec_t(alphabet_size(), 0)),
    lambda_(periodic, 1), depth_(parent == NULL ? 0 : parent->depth_+1),
    count2_added_(false)
{
}

void node::copy(node* N)
{
  // copy normal field
  *this = *N;
  // copy pointer field
  for (int i = 0; i < alphabet_size(); i++) {
    char c = pos2ch(i);
    if (N->children_[i] != NULL) { // if child is not NULL
      children_[i] = new node(alphabet_, periodic(), c, this); 
      children_[i]->copy(N->children_[i]);
      children_[i]->parent_ = this;
    }
  }
}

vec_t node::prob_imm(int frame)
{
  node* N = this;
  double lambda = 1;
  vec_t f = zero_vec_t(alphabet_size());
  while (N != NULL) {
    f += lambda * N->lambda_[frame] * N->f_[frame];
    lambda *= 1-N->lambda_[frame];
    N = N->parent_;
  }
  return f;
}

str_t node::str()
{
  node* N = this;
  str_t s;
  s.reserve(N->depth_);
  while (N != NULL && N->parent_ != NULL) {
    s = s + N->ch_;
    N = N->parent_;
  }
  return s;
}

vom::vom(const str_t& alphabet, int periodic) 
  : alphabet_(alphabet), periodic_(periodic), 
    root_(NULL),
    count_frame_(periodic), depth_frame_(periodic)
{
}


// ==================== Visitors =======================

struct create_visitor {
  int order_;
  create_visitor(int order) : order_(order) {}
  void operator() (node* N) { // N must be not NULL here 
    if (N == NULL) return;
    if (N->depth_ < order_) {
      for (int j = 0; j < N->alphabet_size(); j++) {
	if ( (N->children_[j] = new node(N->alphabet_, N->periodic_, N->pos2ch(j), N) ) == NULL ) {
	  FATAL_ERROR("Cannot allocate memory");
	}	
      }
    }
  }
};

struct addcount_visitor {
  int frame_;
  int chpos_;
  bool continue_;
  addcount_visitor(int frame, int chpos) : frame_(frame), chpos_(chpos), continue_(true) {}
  void operator() (node* N) { // N must be not NULL here 
    if (N == NULL) return;
    N->count_[frame_][chpos_] += 1;    
  }
};

struct addcount1_visitor {
  int frame_;
  int chpos_;
  bool continue_;
  addcount1_visitor(int frame, int chpos) : frame_(frame), chpos_(chpos), continue_(true) {}
  void operator() (node* N) { // N must be not NULL here 
    if (N == NULL) return;
    N->count1_[frame_][chpos_] += 1;    
  }
};

struct normalize_visitor {
  double scale_;
  normalize_visitor(double scale) : scale_(scale) {}
  void operator() (node* N) { // N must be not NULL here 
    if (N == NULL) return;
    for (int frame = 0; frame < N->periodic(); frame++) {
      const vec_t& count = N->count_[frame];
      vec_t& ft = N->f_[frame];
      bool has_zero = false;
      double min = DOUBLE_INFINITY;
      for (int i = 0; i < N->alphabet_size(); i++) {
	if (count(i) == 0) has_zero = true;
	else if (count(i) < min) min = count(i);
      }
      if (min == DOUBLE_INFINITY) min = 1;
      if (!has_zero) 
	ft = count;
      else 
	ft = count + scalar_vec_t(N->alphabet_size(), min*scale_);
      ft /= ublas::sum(ft);
    }
  }
};

struct probability_visitor {
  int frame_;
  int chpos_;
  bool continue_;
  double result_;
  probability_visitor(int frame, int chpos) 
    : frame_(frame), chpos_(chpos), continue_(true), result_(0) {}
  void operator() (const node* N) { // N must be not NULL here 
    if (N == NULL || !N->has_table(frame_)) {
      continue_ = false;
      return;
    }
    result_ = N->f_[frame_][chpos_];
  }
};

struct probability_lambda_visitor {
  int frame_;
  int chpos_;
  bool continue_;
  double result_;
  probability_lambda_visitor(int frame, int chpos) 
    : frame_(frame), chpos_(chpos), continue_(true), result_(0) {}
  void operator() (const node* N) { // N must be not NULL here 
    double lambda = N->lambda_[frame_];
    if (N == NULL || !N->has_table(frame_) || lambda == 0) {
      continue_ = false;
      return;
    }
    result_ = lambda * N->f_[frame_][chpos_] + (1-lambda) * result_;
  }
};

struct statistics_visitor {
  vom& model_;
  statistics_visitor(vom& model) : model_(model) {}
  void operator() (node* N) { 
    if (N == NULL) return;
    model_.count_ ++;
    if (N->depth_ > model_.depth_) model_.depth_ = N->depth_;
    for (int frame = 0; frame < model_.periodic_; frame++) {
      if (N->has_table(frame)) {
	model_.count_frame_[frame]++;
	if (N->depth_ > model_.depth_frame_[frame]) model_.depth_frame_[frame] = N->depth_;
      }
    }
  }
};

struct prune_visitor {
  double eta_;

  prune_visitor(double eta) : eta_(eta) {}
  void operator() (node* N) { 
    if (N == NULL) return;
    for (int frame = 0; frame < N->periodic_; frame++) {
      bool can_prune = N->parent_ != NULL && N->parent_->has_table(frame) 
	&& N->has_table(frame);
      BOOST_FOREACH(node* child, N->children_) 
	if (child && child->has_table(frame)) can_prune = false;
      if (!can_prune) continue;
      
      double KL = CommonFunction::KLdivergence(N->f_[frame], N->parent_->f_[frame]);
      if (KL < eta_) N->f_[frame] = vec_t();
    }
  }
};

struct prune_lambda_visitor {
  typedef std::vector<vec_t> data_t;
  double eta_;

  prune_lambda_visitor(double eta) : eta_(eta) {}
  void operator() (node* N, const data_t& data) { 
    if (N == NULL) return;
    for (int frame = 0; frame < N->periodic_; frame++) {
      bool can_prune = N->parent_ != NULL && N->parent_->has_table(frame) && N->has_table(frame);
      BOOST_FOREACH(node* child, N->children_) 
	if (child && child->has_table(frame) && child->lambda_[frame] > 0) can_prune = false;
      if (!can_prune) continue;
      
      const double& lambda = N->lambda_[frame];
      const vec_t& p_imm = data[frame];
      const vec_t& p = N->f_[frame];
      double KL = 0; //CommonFunction::KLdivergence((*N)->f_[pr], parent->f_[pr]);
      for (int i = 0; i < N->alphabet_size(); i++)
	KL += p(i) * log(lambda * p(i) / p_imm(i) + 1-lambda);

      if (KL < eta_ || lambda == 0) {
	N->f_[frame] = vec_t();
	N->lambda_[frame] = 0;
      }
    }
  }

  data_t compute_data(node* N, const data_t& data) {
    data_t new_data = data;
    for (int frame = 0; frame < N->periodic_; frame++) {
      const double& lambda = N->lambda_[frame];
      if (N->has_table(frame)) {
	const vec_t& ft = N->f_[frame];
	new_data[frame] = lambda*ft + (1-lambda)*data[frame];
      }
    }
    return new_data;
  }
};

struct chisquared_visitor {
  typedef std::vector<vec_t> data_t;
  int T_;
  boost::math::chi_squared chi_dist;
  chisquared_visitor(int T) : T_(T), chi_dist(3) {}
  data_t operator() (node* N, const data_t& data) { 
    if (N == NULL) return data;
    data_t new_data = data;

    for (int frame = 0; frame < N->periodic_; frame++) {
      double& lambda = N->lambda_[frame];
      const vec_t& imm_ft = data[frame];
      const vec_t& count = N->count_[frame];
      const vec_t& ft = N->f_[frame];
      double total_count = ublas::sum(count);

      if (!N->has_table(frame)) {
	lambda = 0;
	continue;
      }
      else if (N->parent_ == NULL || total_count > T_) 
	lambda = 1;
      else {
	double xi = 0;
	for (int j = 0; j < N->alphabet_size(); j++) 
	  xi += ((count(j) == 0) ? (imm_ft[j]*total_count) : ((sqr(count(j) - imm_ft[j]*total_count)) / (imm_ft[j]*total_count)));
	double q = boost::math::cdf(chi_dist, xi);
	lambda = (q < 0.5) ? 0 : q * total_count / T_;
      }

      new_data[frame] = lambda*ft + (1-lambda)*data[frame];
    }

    return new_data;
  }
};

struct statistics_lambda_visitor {
  vom& model_;
  statistics_lambda_visitor(vom& model) : model_(model) {}
  void operator() (node* N) { 
    if (N == NULL) return;
    model_.count_ ++;
    if (N->depth_ > model_.depth_) model_.depth_ = N->depth_;
    for (int frame = 0; frame < model_.periodic_; frame++) {
      if (N->has_table(frame) && N->lambda_[frame] > 0) {
	model_.count_frame_[frame]++;
	if (N->depth_ > model_.depth_frame_[frame]) model_.depth_frame_[frame] = N->depth_;
      }
    }
  }
};

// ===================== Traverse Algorithms ========================
struct TreeTraverse {
  template <typename VISITOR>
  static void dfs(node* N, VISITOR& v) {
    if (N == NULL) return; // reach leaf
    v(N);
    BOOST_FOREACH(node* child, N->children_) dfs(child, v);
  }

  template <typename VISITOR, typename DATA>
  static void dfs(node* N, VISITOR& v, const DATA& data) {
    if (N == NULL) return; // reach leaf
    DATA new_data = v(N, data);
    BOOST_FOREACH(node* child, N->children_) dfs(child, v, new_data);
  }

  template <typename VISITOR>
  static void dfs_backorder(node* N, VISITOR& v) {
    if (N == NULL) return; // reach leaf
    BOOST_FOREACH(node* child, N->children_) dfs_backorder(child, v);
    v(N);
  }

  template <typename VISITOR, typename DATA>
  static void dfs_backorder(node* N, VISITOR& v, const DATA& data) {
    if (N == NULL) return; // reach leaf
    DATA new_data = v.compute_data(N, data);
    BOOST_FOREACH(node* child, N->children_) dfs_backorder(child, v, new_data);
    v(N, data);
  }

  template <typename VISITOR>
  static void reach_down(node* N, const str_t& s, int p, VISITOR& v) {
    if (N == NULL) return;
    v(N);
    if (!v.continue_) return; // terminate signaled
    p--;
    if (p < 0) return;
    reach_down(N->children_[N->ch2pos(s[p])], s, p, v);
  }
};

// ========================= VOM ====================================

void vom::build_fom_flat(const Sequences& train, int order)
{
  root_ = new node(&alphabet_, periodic_, 'X', NULL);

  // create tree
  create_visitor create(order);
  TreeTraverse::dfs(root_, create);

  BOOST_FOREACH(const str_t& s, train) {
    for (int i = 0; i < sz(s); i++) {
      addcount_visitor add_count(i % periodic_, root_->ch2pos(s[i]));
      TreeTraverse::reach_down(root_, s, i, add_count);
    }
  }

  normalize_visitor normalize(0.1);
  TreeTraverse::dfs(root_, normalize);
}

struct addcount2_visitor {
  int frame_;
  int chpos_;
  bool continue_;
  addcount2_visitor(int frame, int chpos) : frame_(frame), chpos_(chpos), continue_(true) {}
  void operator() (node* N) { // N must be not NULL here 
    if (N == NULL || N->count2_added_) return; // only add count at leaves
    N->count_[frame_][chpos_] += 1;
  }
};

struct checkcount_visitor {
  int frame_;
  bool continue_;
  int level_;
  double count_;
  checkcount_visitor(int frame, int level) : frame_(frame), continue_(true), level_(level), count_(-1) {}
  void operator() (node* N) { // N must be not NULL here 
    if (N == NULL) return; 
    if (N->depth_ == level_) count_ = ublas::sum(N->count_[frame_]);
  }
};

struct normalize_cut_visitor {
  double scale_;
  int Tmin_;
  normalize_cut_visitor(double scale, int Tmin) : scale_(scale), Tmin_(Tmin) {}
  void operator() (node* N) { // N must be not NULL here 
    if (N == NULL) return;
    for (int frame = 0; frame < N->periodic(); frame++) {
      const vec_t& count = N->count_[frame];
      vec_t& ft = N->f_[frame];

      if (ublas::sum(count) < Tmin_) { // cut if not exceed threshold
	ft = vec_t();
	continue;
      }

      bool has_zero = false;
      double min = DOUBLE_INFINITY;
      for (int i = 0; i < N->alphabet_size(); i++) {
	if (count(i) == 0) has_zero = true;
	else if (count(i) < min) min = count(i);
      }
      if (min == DOUBLE_INFINITY) min = 1;
      if (!has_zero) 
	ft = count;
      else 
	ft = count + scalar_vec_t(N->alphabet_size(), min*scale_);
      ft /= ublas::sum(ft);
    }
  }
};

struct node_str_comp {
  bool operator() (node* N, const str_t& s) const {
    str_t N_s = N->str();
    for (int i = sz(s)-1; i >= 0; i--)
      if (N_s[i] < s[i]) return true;
      else if (N_s[i] > s[i]) return false;
    return false;
  }
  
  bool equal (node* N, const str_t& s) const {
    str_t N_s = N->str();
    return N_s.compare(s) == 0;
  }
};

void vom::build_vom_flat(const Sequences& train, int Tmin) 
{
  root_ = new node(&alphabet_, periodic_, 'X', NULL);
  node_vec_t lv_, next_lv_;
  node_vec_t *lv = &lv_, *next_lv = &next_lv_;

  lv->pb(root_);
  int count_lv = 0;

  // prepare idx
  std::vector<vi_t> idx_(sz(train)), next_idx_(sz(train));
  std::vector<vi_t> *idx = &idx_, *next_idx = &next_idx_;
  int sum_idx = 0;
  for (int s_idx = 0; s_idx < sz(train); s_idx++) {
    (*idx)[s_idx].reserve(sz(train[s_idx]));
    for (int i = 0; i < sz(train[s_idx]); i++) (*idx)[s_idx].pb(i);
    sum_idx += sz(train[s_idx]);
  }

  while (sz(*lv) > 0) {
    MSG("lv=" << (count_lv) << " sz(lv)=" << sz(*lv)  << " sum_idx=" << sum_idx);

    /*
    // this is brute force
    BOOST_FOREACH(const str_t& s, train) {
      for (int i = 0; i < sz(s); i++) {
	addcount2_visitor add_count2(i % periodic_, root_->ch2pos(s[i]));
	TreeTraverse::reach_down(root_, s, i, add_count2);
      }
    }
    BOOST_FOREACH(node* N, *lv) N->count2_added_ = true;
    */

    // add count
    for (int s_idx = 0; s_idx < sz(train); s_idx++) {
      const str_t& s = train[s_idx];
      int current_size = sz((*idx)[s_idx]);
      for (int i_idx = 0; i_idx < current_size; i_idx++) {
	int i = (*idx)[s_idx][i_idx];
	addcount2_visitor add_count2(i % periodic_, root_->ch2pos(s[i]));
	TreeTraverse::reach_down(root_, s, i, add_count2);
      }
    }
    BOOST_FOREACH(node* N, *lv) N->count2_added_ = true;
    // prepare next idx
    sum_idx = 0;
    for (int s_idx = 0; s_idx < sz(train); s_idx++) {
      const str_t& s = train[s_idx];
      int current_size = sz((*idx)[s_idx]);
      (*next_idx)[s_idx].clear();
      (*next_idx)[s_idx].reserve(current_size);
      for (int i_idx = 0; i_idx < current_size; i_idx++) {
	int i = (*idx)[s_idx][i_idx];
	checkcount_visitor check_count(i % periodic_, count_lv);
	TreeTraverse::reach_down(root_, s, i, check_count);
	if (check_count.count_ >= Tmin) (*next_idx)[s_idx].pb(i);
      }
      sum_idx += sz((*next_idx)[s_idx]);
    }
    std::swap(idx, next_idx);

    /* 
    // a more clever way, utilize the fact that nodes in each level are sorted by their context
    BOOST_FOREACH(const str_t& s, train) {
      if (count_lv == 0) {       // root case
	for (int i = 0; i < sz(s); i++) {
	  int frame = i % periodic_;
	  int chpos = root_->ch2pos(s[i]);
	  root_->count_[frame][chpos]++;
	}
      }
      else { 
	node_str_comp comp;
	typedef node_vec_t::iterator iterator;
	for (int i = count_lv; i < sz(s); i++) { // bisection to find the node
	  int frame = i % periodic_;
	  int chpos = root_->ch2pos(s[i]);
	  str_t sub_s = s.substr(i-count_lv, count_lv);
	  iterator it = std::lower_bound(all_element(*lv), sub_s, comp); 
	  if (it != lv->end() && comp.equal(*it, sub_s)) {
	    node* N = *it;
	    N->count_[frame][chpos]++;
	  }
	}
      }
    }
    */

    // populate next level
    next_lv->clear();
    next_lv->reserve(sz(*lv)*sz(alphabet_));
    BOOST_FOREACH(node* N, *lv) {
      bool exceed_threshold = false;
      for (int frame = 0; frame < periodic_; frame++) 
	if (ublas::sum(N->count_[frame]) >= Tmin) exceed_threshold = true;
      
      if (exceed_threshold) {
	for (int j = 0; j < sz(alphabet_); j++) {
	  N->children_[j] = new node(N->alphabet_, N->periodic_, N->pos2ch(j), N);
	  next_lv->pb(N->children_[j]);
	}
      }
    }
    std::swap(lv, next_lv);
    count_lv++;
  }

  normalize_cut_visitor normalize(0.1, Tmin);
  TreeTraverse::dfs(root_, normalize);  
}

double vom::log_likelihood(const str_t& s)
{
  double ll = 0;
  for (int i = 0; i < sz(s); i++) {
    probability_visitor probability(i % periodic_, root_->ch2pos(s[i]));
    TreeTraverse::reach_down(root_, s, i, probability);
    ll += log(probability.result_);
  }
  return ll;
}

double vom::log_likelihood_lambda(const str_t& s)
{
  double ll = 0;
  for (int i = 0; i < sz(s); i++) {
    probability_lambda_visitor probability(i % periodic_, root_->ch2pos(s[i]));
    TreeTraverse::reach_down(root_, s, i, probability);
    ll += log(probability.result_);
  }
  return ll;
}

void vom::compute_stats()
{
  count_ = 0; 
  depth_ = 0;
  std::fill(all_element(count_frame_), 0);
  std::fill(all_element(depth_frame_), 0);
  statistics_visitor statistics(*this);
  TreeTraverse::dfs(root_, statistics);
}

void vom::print_stats(const str_t& outFile)
{
  FileAccessWrite(out, outFile);
  out << " " << std::setw(10) << count_ 
      << " " << std::setw(10) << depth_;
  for (int i = 0; i < periodic_; i++)
    out << " " << std::setw(10) << count_frame_[i] 
	<< " " << std::setw(10) << depth_frame_[i]; 
  out << ENDL;
}

void vom::prune(double eta) 
{
  prune_visitor prune_v(eta);
  TreeTraverse::dfs_backorder(root_, prune_v);
}

void vom::prune_lambda(double eta) 
{
  std::vector<vec_t> imm_f(periodic_, zero_vec_t(sz(alphabet_)));
  prune_lambda_visitor prune_v(eta);
  MSG("0");
  TreeTraverse::dfs_backorder(root_, prune_v, imm_f);
  MSG("1");
}

void vom::copy(vom& m)
{
  // copy normal field
  *this = m;
  // pointer field
  root_ = new node(&alphabet_, periodic_, 'X', NULL);
  root_->copy(m.root_);
}

// ==================== Interpolated Markov Models (chi-squared) ================
void vom::compute_lambda_chi_square(int T)
{
  std::vector<vec_t> imm_f(periodic_, zero_vec_t(sz(alphabet_)));
  chisquared_visitor chisquare(T);
  TreeTraverse::dfs(root_, chisquare, imm_f);
}

void vom::compute_stats_lambda()
{
  count_ = 0; 
  depth_ = 0;
  std::fill(all_element(count_frame_), 0);
  std::fill(all_element(depth_frame_), 0);
  statistics_lambda_visitor statistics(*this);
  TreeTraverse::dfs(root_, statistics);
}

// ==================  Deleted Interpolation ========================
void vom::populate_next_level(node_vec_t* lv, node_vec_t* next_lv) {
  next_lv->clear();
  next_lv->reserve(sz(*lv) * sz(alphabet_));
  BOOST_FOREACH(node* N, *lv) {
    BOOST_FOREACH(node* child, N->children_) if (child != NULL) next_lv->pb(child);
  }
}

struct compare_node {
  int frame_;
  bool min_;
  compare_node(int frame, bool min) : frame_(frame), min_(min) {}
  bool operator() (node* N1, node* N2) const {
    double n1 = ublas::sum(N1->count_[frame_]), 
      n2 = ublas::sum(N2->count_[frame_]);
    if (min_) 
      return n1 < n2 && n1 > 1e-10;
    else return n1 > n2;
  }
};

void vom::split_to_bin(node_vec_t* lv, std::vector<node_vec_t>& bins, int frame, double delta) {
  typedef node_vec_t::iterator iterator;

  compare_node min_comp(frame, true), max_comp(frame, false);
  iterator min = std::min_element(all_element(*lv), min_comp),
    max = std::min_element(all_element(*lv), max_comp);
  double minN = ublas::sum((*min)->count_[frame]), 
    maxN = ublas::sum((*max)->count_[frame]);

  bins.clear();
  if (maxN == 0) return;
  if (minN == 0) return;

  int M = (int) (floor( log(maxN/minN) / log(delta) ) + 1);
  bins.resize(M);
  BOOST_FOREACH(node* N, *lv) {
    double n = ublas::sum(N->count_[frame]);
    if (n == 0) N->lambda_[frame] = 0;
    else {
      int bin_idx = (int) (floor( log(n/minN) / log(delta) ));
      bins[bin_idx].pb(N);
    }
  }

  MSG("M=" << M);
  for (int i = 0; i < M; i++) std::cerr << std::setw(8) << sz(bins[i]);
  std::cerr << ENDL;
}

double vom::bin_func(node_vec_t& bin, int frame, double lambda)
{
  double s = 0;
  BOOST_FOREACH(node* N, bin) {
    if (!N->has_table(frame)) continue;
    node* parent = N->parent_;
    vec_t p_imm = parent->prob_imm(frame);
    const vec_t& ft = N->f_[frame];
    vec_t diff = ft-p_imm;
    const vec_t& count = N->count1_[frame];
    //    std::cerr << "ch=" << N->ch_ << ENDL;
    //    std::cerr << "count=" << count << ENDL << " diff=" << diff << ENDL << "p_imm=" << p_imm << ENDL;
    for (int j = 0; j < N->alphabet_size(); j++)
      s += count(j) * diff(j) / (lambda * diff(j) + p_imm(j));
  }
  return s;
}

double vom::optimize_bin(node_vec_t& bin, int frame)
{
  double f0 = bin_func(bin, frame, 0);
  double f1 = bin_func(bin, frame, 1);
  //  std::cerr << std::setw(8) << std::fixed << std::setprecision(2) << f0;
  //  std::cerr << std::setw(8) << std::fixed << std::setprecision(2) << f1;

  if (f0 <= 0) return 0; // special case, empty bin return 0
  if (f1 >= 0) return 1;

  // bisection, start f(0) > 0, f(1) < 0
  double l = 0, r = 1;
  while (fabs(r-l) > 1e-8) {
    double m = (l+r)/2;
    //MSG("bi " << l << " " << r << " " << m);
    double f_m = bin_func(bin, frame, m);
    if (f_m <= 0) r = m;
    else l = m;    
  }
  return l;
}

void vom::compute_lambda_deleted_interpolation(const Sequences& train, double delta)
{
  BOOST_FOREACH(const str_t& s, train) {
    for (int i = 0; i < sz(s); i++) {
      addcount1_visitor add_count(i % periodic_, root_->ch2pos(s[i]));
      TreeTraverse::reach_down(root_, s, i, add_count);
    }
  }

  node_vec_t lv_, next_lv_;
  node_vec_t *lv = &lv_, *next_lv = &next_lv_;
  
  std::fill(all_element(root_->lambda_), 1.0);

  lv->pb(root_); populate_next_level(lv, next_lv);
  std::swap(lv, next_lv);

  int n_lv = 1;
  while (sz(*lv) > 0) { // process each level
    MSG("lv=" << (n_lv++) << " size(lv)=" << sz(*lv));
    for (int frame = 0; frame < periodic_; frame++) {
      std::vector<node_vec_t> bins;
      split_to_bin(lv, bins, frame, delta);
      
      BOOST_FOREACH(node_vec_t& bin, bins) {
	double lambda = optimize_bin(bin, frame);
	std::cerr << std::setw(8) << std::fixed << std::setprecision(2) << lambda;
	BOOST_FOREACH(node* N, bin) N->lambda_[frame] = lambda;
      }
      std::cerr << ENDL;
    }

    populate_next_level(lv, next_lv);
    std::swap(lv, next_lv);
  }
}



VOM_END
